from math import *
import sympy

def eqReid78(lambda0,xi,D):
    beta = 1 / sqrt(2) * sqrt(1-1/lambda0 *(1 - 2 / (xi**2 * D**2)))
    varphi_1 = acos(1/(beta * xi * D * sqrt(lambda0)))
    c_D = (2 * sympy.elliptic_e(varphi_1,beta**2) \
           - sympy.elliptic_f(varphi_1,beta**2) ) \
            / (D * xi * sqrt(lambda0)) 
    delta_D = sqrt(1 - 1 / lambda0**2) - 2*c_D
    return delta_D

# 利用eqReid78()输出一系列力-位移
def from_F_get_U(U_end,step,Y,E_p,D,t):
    xi = sqrt(12.0 * Y / (D * E_p * t))
    U = [0]
    F = [1]
    while U[-1] < U_end:
        F0 = F[-1] + step
        F.append(F0)
        U.append(eqReid78(F0,xi,D))
    return (U,F)

# 返还量纲
def get_trueUF(Uno,Fno,P_0,D):
    trueF = [F0 * P_0 for F0 in Fno]
    trueU = [U0 * D for U0 in Uno]
    return (trueU,trueF)   

#将结果写入文件
def write2file(path,Col1,Col2):
    f = open(path,'w')
    n = len(Col1)
    for i in range(0,n):
        line = str(Col1[i]) + '    ' + str(Col2[i]) + '\n'
        f.write(line)
    f.close()

if __name__ == "__main__":
    Y = 160 # 屈服应力
    D = 130 # 直径
    E_p = 758 # 硬化模量
    L = 40 # 管长
    t = 3 # 壁厚
    U_end = 0.8 # 压溃效率，理论取值范围（0,0.98），实际上0.9以上计算就很吃力了
    P_0 = 2 * Y * t**2 * L / D # 初始塌陷荷载
    (U,F) = from_F_get_U(U_end,0.01,Y,E_p,D,t)
    (tU,tF) = get_trueUF(U,F,P_0,D)
    print(len(tF),len(tF)==len(tU)) # 输出数据量并核对两列数据的长度是否相等
    path = ".\\test.txt" # 输出文件的路径
    write2file(path,tU,tF) # 写入数据